function [Uout,emax] = propag2Q(Uop,t,nch,np,fieldeps,emax,E0,muop)
       
%        propagation subroutine by Chebychev reccursion

%        calculating the spectral range of the Hamiltonian operator
        emax = sqrt((max(E0)-min(E0))^2+fieldeps^2);
        
        range=t*emax/2.d0;

%        phase correction
        Phase = exp(i*range);


%        prepairing the exapnsion coefficients

       bj=   besselj(0:nch-1,range);

        
        Uop1 = Uop;

%        the first chebychev polynomial J(0)

        Uop = Uop*bj(1)*Phase;

%        The second Chebychev polynomial

        Uop2 =  multi2Q(Uop1,np,emax,fieldeps,E0,muop);

        Uop = Uop + Uop2*bj(2)*2.d0*i*Phase;

%        The Chebychev recuursion and accumulation of results:

        for jj = 2:nch-1

%        ***   compute cii=(0,1)**j ***
        j4 = mod(jj,4);
        cii = i^j4;
        Uop3 =  multi2Q(Uop2,np,emax,fieldeps,E0,muop) ;

%        Chebychev reccursion relation

        Uop3 = 2.d0*Uop3 - Uop1;
 
        Uop1 = Uop2;
        Uop2 = Uop3;

%        accumulation of result

        Uop = Uop + Uop2*bj(jj+1)*2.d0*cii*Phase;
        end
    Uout = Uop;
        end 